﻿/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://github.com/shenyczz/CSharpKit
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime;
using System.Text;
using CSharpKit.Numerics.Theories;
using CSharpKit.Threading.Tasks;

namespace CSharpKit.Numerics.LinearAlgebra
{
    /// <summary>
    /// 矩阵 <para/>
    /// </summary>
    /// <typeparam name="T"></typeparam>
    /// <remarks>
    /// 在矩阵中，若数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵；
    /// 与之相反，若非0元素数目占大多数时，则称该矩阵为稠密矩阵。定义非零元素的总数比上矩阵所有元素的总数为矩阵的稠密度。<br/>
    /// 方阵 - square matrix <br/>
    /// 正交矩阵 - orthogonal matrix (columns are orthogonal unit vectors meaning QTQ = I)<br/>
    /// 上三角阵 - upper triangular matrix (also called right triangular matrix R)
    /// 下三角阵 - lower triangular matrix (also called left  triangular matrix L)
    /// </remarks>
    [Serializable]
    public abstract partial class Matrix<T> : Numeric<T>, IFormattable, ICloneable<Matrix<T>>, IEquatable<Matrix<T>>
        where T : struct, IFormattable, IEquatable<T>
    {
        #region Constructors

        protected Matrix(MatrixStorage<T> storage)
        {
            Storage = storage;
            RowCount = storage.RowCount;
            ColumnCount = storage.ColumnCount;
        }

        #endregion

        /// <summary>
        /// 存储序（默认列主序）
        /// </summary>
        /// <value></value>
        public static MajorOrders MajorOrder { get; set; } = MajorOrders.Column;

        public static MatrixBuilder<T> Builder => BuilderInstance<T>.MatrixBuilder;


        #region Properties

        public int ColumnCount { get; private set; }

        public int RowCount { get; private set; }

        public MatrixStorage<T> Storage { get; private set; }

        #endregion


        #region Ising

        // 矩阵特征属性

        // 方块矩阵 Square Matrices
        // 对称矩阵 Symmetric Matrices
        // 单位矩阵 Identity Matrices
        // 艾米矩阵 Hermitian Matrices

        /// <summary>
        /// 是否对称矩阵
        /// 对称矩阵（Symmetric Matrices）是指以主对角线为对称轴，各元素对应相等的矩阵。
        /// 在线性代数中，对称矩阵是一个方形矩阵，其转置矩阵和自身相等。
        /// </summary>
        /// <returns></returns>
        public virtual bool IsSymmetric()
        {
            if (RowCount != ColumnCount)
            {
                return false;
            }

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = row + 1; column < ColumnCount; column++)
                {
                    if (!At(row, column).Equals(At(column, row)))
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        public virtual bool IsHermitian() => throw new Exception();

        #endregion


        #region Indexer - 索引器

        public T this[int row, int column]
        {
            get { return At(row, column); }
            set { At(row, column, value); }
        }

        [TargetedPatchingOptOut("Performance critical to inline across NGen image boundaries")]
        internal T At(int row, int column)
        {
            return Storage.At(row, column);
        }
        [TargetedPatchingOptOut("Performance critical to inline across NGen image boundaries")]
        internal void At(int row, int column, T value)
        {
            Storage.At(row, column, value);
        }

        #endregion


        #region ICloneable<Matrix<T>>

        public Matrix<T> Clone()
        {
            var result = Builder.SameAs(this);
            Storage.CopyToUnchecked(result.Storage, ExistingData.AssumeZeros);
            return result;
        }

        #endregion

        #region IEquatable<Matrix<T>>

        public bool Equals(Matrix<T> other)
        {
            return other != null
                && Storage.Equals(other.Storage);
        }

        public override bool Equals(object obj)
        {
            return Equals(obj as Matrix<T>);
        }

        public override int GetHashCode()
        {
            return base.GetHashCode()
                ^ Storage.GetHashCode()
                ;
        }

        #endregion

        #region IFormattable

        public sealed override string ToString()
        {
            return ToString("G6", null);
            //return string.Concat(ToTypeString(), Environment.NewLine, ToMatrixString("G6"));
        }

        public string ToString(string format, IFormatProvider formatProvider = null)
        {
            return string.Concat(ToTypeString(), Environment.NewLine, ToMatrixString(format, formatProvider));
        }

        public virtual string ToTypeString()
        {
            return string.Format("{0} {1}x{2}-{3}", GetType().Name, RowCount, ColumnCount, typeof(T).Name);
        }

        public string ToMatrixString(string format, IFormatProvider provider = null)
        {
            if (string.IsNullOrEmpty(format))
            {
                format = "G6";
            }

            return ToMatrixString(8, 4, 5, 2, 76,
                "..", "..", "..",
                "  ", Environment.NewLine,
                x => x.ToString(format, provider));
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="upperRows">上部行数</param>
        /// <param name="lowerRows">下部行数</param>
        /// <param name="minLeftColumns">最小左边列数</param>
        /// <param name="rightColumns">右边列数</param>
        /// <param name="maxWidth">最大宽度</param>
        /// <param name="horizontalEllipsis">水平省略符</param>
        /// <param name="verticalEllipsis">垂直省略符</param>
        /// <param name="diagonalEllipsis">对角线省略符</param>
        /// <param name="columnSeparator">列分隔符</param>
        /// <param name="rowSeparator">行分隔符</param>
        /// <param name="formatValue"></param>
        /// <returns></returns>
        public string ToMatrixString(int upperRows, int lowerRows, int minLeftColumns, int rightColumns, int maxWidth,
            string horizontalEllipsis, string verticalEllipsis, string diagonalEllipsis,
            string columnSeparator, string rowSeparator, Func<T, string> formatValue)
        {
            return FormatStringArrayToString(
                ToMatrixStringArray(upperRows, lowerRows, minLeftColumns, rightColumns, maxWidth, columnSeparator.Length, horizontalEllipsis, verticalEllipsis, diagonalEllipsis, formatValue),
                columnSeparator, rowSeparator);
        }

        /// <summary>
        /// Returns a string 2D array that summarizes the content of this matrix.
        /// </summary>
        public string[,] ToMatrixStringArray(int upperRows, int lowerRows, int minLeftColumns, int rightColumns, int maxWidth, int padding,
            string horizontalEllipsis, string verticalEllipsis, string diagonalEllipsis, Func<T, string> formatValue)
        {
            upperRows = System.Math.Max(upperRows, 1);
            lowerRows = System.Math.Max(lowerRows, 0);
            minLeftColumns = System.Math.Max(minLeftColumns, 1);
            maxWidth = System.Math.Max(maxWidth, 12);

            int upper = RowCount <= upperRows ? RowCount : upperRows;
            int lower = RowCount <= upperRows ? 0 : RowCount <= upperRows + lowerRows ? RowCount - upperRows : lowerRows;
            bool rowEllipsis = RowCount > upper + lower;
            int rows = rowEllipsis ? upper + lower + 1 : upper + lower;

            int left = ColumnCount <= minLeftColumns ? ColumnCount : minLeftColumns;
            int right = ColumnCount <= minLeftColumns ? 0 : ColumnCount <= minLeftColumns + rightColumns ? ColumnCount - minLeftColumns : rightColumns;

            var columnsLeft = new List<Tuple<int, string[]>>();
            for (int j = 0; j < left; j++)
            {
                columnsLeft.Add(FormatColumn(j, rows, upper, lower, rowEllipsis, verticalEllipsis, formatValue));
            }

            var columnsRight = new List<Tuple<int, string[]>>();
            for (int j = 0; j < right; j++)
            {
                columnsRight.Add(FormatColumn(ColumnCount - right + j, rows, upper, lower, rowEllipsis, verticalEllipsis, formatValue));
            }

            int chars = columnsLeft.Sum(t => t.Item1 + padding) + columnsRight.Sum(t => t.Item1 + padding);
            for (int j = left; j < ColumnCount - right; j++)
            {
                var candidate = FormatColumn(j, rows, upper, lower, rowEllipsis, verticalEllipsis, formatValue);
                chars += candidate.Item1 + padding;
                if (chars > maxWidth)
                {
                    break;
                }
                columnsLeft.Add(candidate);
            }

            int cols = columnsLeft.Count + columnsRight.Count;
            bool colEllipsis = ColumnCount > cols;
            if (colEllipsis)
            {
                cols++;
            }

            var array = new string[rows, cols];
            int colIndex = 0;
            foreach (var column in columnsLeft)
            {
                for (int i = 0; i < column.Item2.Length; i++)
                {
                    array[i, colIndex] = column.Item2[i];
                }
                colIndex++;
            }
            if (colEllipsis)
            {
                int rowIndex = 0;
                for (var row = 0; row < upper; row++)
                {
                    array[rowIndex++, colIndex] = horizontalEllipsis;
                }
                if (rowEllipsis)
                {
                    array[rowIndex++, colIndex] = diagonalEllipsis;
                }
                for (var row = RowCount - lower; row < RowCount; row++)
                {
                    array[rowIndex++, colIndex] = horizontalEllipsis;
                }
                colIndex++;
            }
            foreach (var column in columnsRight)
            {
                for (int i = 0; i < column.Item2.Length; i++)
                {
                    array[i, colIndex] = column.Item2[i];
                }
                colIndex++;
            }
            return array;
        }

        /// <summary>
        /// 格式化字符串数组为字符串
        /// </summary>
        /// <param name="array">字符串数组</param>
        /// <param name="colSeparator">列分隔符</param>
        /// <param name="rowSeparator">行分割符</param>
        /// <returns></returns>
        private string FormatStringArrayToString(string[,] array, string colSeparator, string rowSeparator)
        {
            var rows = array.GetLength(0);
            var cols = array.GetLength(1);

            var widths = new int[cols];
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    widths[j] = System.Math.Max(widths[j], array[i, j].Length);
                }
            }

            var sb = new StringBuilder();
            for (int i = 0; i < rows; i++)
            {
                sb.Append(array[i, 0].PadLeft(widths[0]));
                for (int j = 1; j < cols; j++)
                {
                    sb.Append(colSeparator);
                    sb.Append(array[i, j].PadLeft(widths[j]));
                }
                sb.Append(rowSeparator);
            }
            return sb.ToString();
        }

        /// <summary>
        /// 格式化列
        /// </summary>
        /// <param name="column"></param>
        /// <param name="height"></param>
        /// <param name="upper"></param>
        /// <param name="lower"></param>
        /// <param name="withEllipsis"></param>
        /// <param name="ellipsis"></param>
        /// <param name="formatValue"></param>
        /// <returns></returns>
        private Tuple<int, string[]> FormatColumn(int column, int height, int upper, int lower, bool withEllipsis, string ellipsis, Func<T, string> formatValue)
        {
            var c = new string[height];
            int index = 0;
            for (var row = 0; row < upper; row++)
            {
                c[index++] = formatValue(At(row, column));
            }
            if (withEllipsis)
            {
                c[index++] = "";
            }
            for (var row = RowCount - lower; row < RowCount; row++)
            {
                c[index++] = formatValue(At(row, column));
            }
            int w = c.Max(x => x.Length);
            if (withEllipsis)
            {
                c[upper] = ellipsis;
            }
            return new Tuple<int, string[]>(w, c);
        }

        #endregion


        #region Clearing

        public void Clear()
        {
            Storage.Clear();
        }
        public void ClearRow(int rowIndex)
        {
            if ((uint)rowIndex >= (uint)RowCount)
            {
                throw new ArgumentOutOfRangeException(nameof(rowIndex));
            }

            Storage.ClearUnchecked(rowIndex, 1, 0, ColumnCount);
        }
        public void ClearColumn(int columnIndex)
        {
            if ((uint)columnIndex >= (uint)ColumnCount)
            {
                throw new ArgumentOutOfRangeException(nameof(columnIndex));
            }

            Storage.ClearUnchecked(0, RowCount, columnIndex, 1);
        }
        public void ClearRows(params int[] rowIndices)
        {
            Storage.ClearRows(rowIndices);
        }
        public void ClearColumns(params int[] columnIndices)
        {
            Storage.ClearColumns(columnIndices);
        }
        public void ClearSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
        {
            Storage.Clear(rowIndex, rowCount, columnIndex, columnCount);
        }


        public void CoerceZero(Func<T, bool> fun)
        {
            MapInplace(x => (fun(x) ? Zero : x), Zeros.AllowSkip);
        }

        public virtual void CoerceZero(T threshold)
        {
            MapInplace(x => NumOperators.LessThan(x, threshold) ? Zero : x, Zeros.AllowSkip);
        }

        #endregion


        #region Copying

        public void CopyTo(Matrix<T> target)
        {
            if (target == null)
            {
                throw new ArgumentNullException(nameof(target));
            }

            Storage.CopyTo(target.Storage);
        }

        #endregion


        #region Enumerating

        public IEnumerable<T> Enumerate()
        {
            return Storage.Enumerate();
        }
        public IEnumerable<T> Enumerate(Zeros zeros = Zeros.Include)
        {
            switch (zeros)
            {
                case Zeros.AllowSkip:
                    return Storage.EnumerateNonZero();

                default:
                    return Storage.Enumerate();
            }
        }
        public IEnumerable<Tuple<int, int, T>> EnumerateIndexed()
        {
            return Storage.EnumerateIndexed();
        }
        public IEnumerable<Tuple<int, int, T>> EnumerateIndexed(Zeros zeros = Zeros.Include)
        {
            switch (zeros)
            {
                case Zeros.AllowSkip:
                    return Storage.EnumerateNonZeroIndexed();
                default:
                    return Storage.EnumerateIndexed();
            }
        }

        public IEnumerable<Vector<T>> EnumerateColumns()
        {
            for (var i = 0; i < ColumnCount; i++)
            {
                yield return Column(i);
            }
        }
        public IEnumerable<Vector<T>> EnumerateColumns(int index, int length)
        {
            var maxIndex = System.Math.Min(index + length, ColumnCount);
            for (var i = System.Math.Max(index, 0); i < maxIndex; i++)
            {
                yield return Column(i);
            }
        }

        public IEnumerable<Tuple<int, Vector<T>>> EnumerateColumnsIndexed()
        {
            for (var i = 0; i < ColumnCount; i++)
            {
                yield return new Tuple<int, Vector<T>>(i, Column(i));
            }
        }
        public IEnumerable<Tuple<int, Vector<T>>> EnumerateColumnsIndexed(int index, int length)
        {
            var maxIndex = System.Math.Min(index + length, ColumnCount);
            for (var i = System.Math.Max(index, 0); i < maxIndex; i++)
            {
                yield return new Tuple<int, Vector<T>>(i, Column(i));
            }
        }

        public IEnumerable<Vector<T>> EnumerateRows()
        {
            for (var i = 0; i < RowCount; i++)
            {
                yield return Row(i);
            }
        }
        public IEnumerable<Vector<T>> EnumerateRows(int index, int length)
        {
            var maxIndex = System.Math.Min(index + length, RowCount);
            for (var i = System.Math.Max(index, 0); i < maxIndex; i++)
            {
                yield return Row(i);
            }
        }
        public IEnumerable<Tuple<int, Vector<T>>> EnumerateRowsIndexed()
        {
            for (var i = 0; i < RowCount; i++)
            {
                yield return new Tuple<int, Vector<T>>(i, Row(i));
            }
        }
        public IEnumerable<Tuple<int, Vector<T>>> EnumerateRowsIndexed(int index, int length)
        {
            var maxIndex = System.Math.Min(index + length, RowCount);
            for (var i = System.Math.Max(index, 0); i < maxIndex; i++)
            {
                yield return new Tuple<int, Vector<T>>(i, Row(i));
            }
        }

        #endregion


        #region Extracting

        public T[,] ToArray()
        {
            return Storage.ToArray();
        }

        public T[][] ToRowArrays()
        {
            return Storage.ToRowArrays();
        }
        public T[][] ToColumnArrays()
        {
            return Storage.ToColumnArrays();
        }

        public T[] ToRowMajorArray()
        {
            return Storage.ToRowMajorArray();
        }
        public T[] ToColumnMajorArray()
        {
            return Storage.ToColumnMajorArray();
        }

        #endregion

        #region Extracting - Vector<T>

        public Vector<T> Row(int index)
        {
            if ((uint)index >= (uint)RowCount)
            {
                throw new ArgumentOutOfRangeException(nameof(index));
            }

            var ret = Vector<T>.Builder.SameAs(this, ColumnCount);
            Storage.CopySubRowToUnchecked(ret.Storage, index, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            return ret;
        }

        public void Row(int index, Vector<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            Storage.CopyRowTo(result.Storage, index);
        }

        public Vector<T> Row(int rowIndex, int columnIndex, int length)
        {
            var ret = Vector<T>.Builder.SameAs(this, length);
            Storage.CopySubRowTo(ret.Storage, rowIndex, columnIndex, 0, length);
            return ret;
        }

        public void Row(int rowIndex, int columnIndex, int length, Vector<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            Storage.CopySubRowTo(result.Storage, rowIndex, columnIndex, 0, length);
        }


        public Vector<T> Column(int index)
        {
            if ((uint)index >= (uint)ColumnCount)
            {
                throw new ArgumentOutOfRangeException(nameof(index));
            }

            var ret = Vector<T>.Builder.SameAs(this, RowCount);
            Storage.CopySubColumnToUnchecked(ret.Storage, index, 0, 0, RowCount, ExistingData.AssumeZeros);
            return ret;
        }

        public void Column(int index, Vector<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            Storage.CopyColumnTo(result.Storage, index);
        }

        public Vector<T> Column(int columnIndex, int rowIndex, int length)
        {
            var ret = Vector<T>.Builder.SameAs(this, length);
            Storage.CopySubColumnTo(ret.Storage, columnIndex, rowIndex, 0, length);
            return ret;
        }

        public void Column(int columnIndex, int rowIndex, int length, Vector<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            Storage.CopySubColumnTo(result.Storage, columnIndex, rowIndex, 0, length);
        }


        public virtual Vector<T> Diagonal()
        {
            var min = System.Math.Min(RowCount, ColumnCount);
            var diagonal = Vector<T>.Builder.SameAs(this, min);

            for (var i = 0; i < min; i++)
            {
                diagonal.At(i, At(i, i));
            }

            return diagonal;
        }


        #endregion

        #region Extracting - Matrix<T>

        /// <summary>
        /// 下三角阵
        /// </summary>
        /// <returns></returns>
        public virtual Matrix<T> LowerTriangle()
        {
            var result = Builder.SameAs(this);

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = 0; column <= row && column < ColumnCount; column++)
                {
                    result.At(row, column, At(row, column));
                }
            }

            return result;
        }
        public virtual void LowerTriangle(Matrix<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = 0; column < ColumnCount; column++)
                {
                    result.At(row, column, row >= column ? At(row, column) : Zero);
                }
            }
        }

        /// <summary>
        /// 上三角阵
        /// </summary>
        public virtual Matrix<T> UpperTriangle()
        {
            var result = Builder.SameAs(this);

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = row; column < ColumnCount; column++)
                {
                    result.At(row, column, At(row, column));
                }
            }

            return result;
        }
        public virtual void UpperTriangle(Matrix<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = 0; column < ColumnCount; column++)
                {
                    result.At(row, column, row <= column ? At(row, column) : Zero);
                }
            }
        }

        public virtual Matrix<T> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
        {
            var result = Builder.SameAs(this, rowCount, columnCount);
            Storage.CopySubMatrixTo(result.Storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, ExistingData.AssumeZeros);
            return result;
        }

        /// <summary>
        /// 严格的下三角形（不包含对角线上元素）
        /// </summary>
        /// <returns>The lower triangle of this matrix.</returns>
        public virtual Matrix<T> StrictlyLowerTriangle()
        {
            var result = Builder.SameAs(this);
            for (var row = 0; row < RowCount; row++)
            {
                var columns = System.Math.Min(row, ColumnCount);
                for (var column = 0; column < columns; column++)
                {
                    result.At(row, column, At(row, column));
                }
            }
            return result;
        }
        public virtual void StrictlyLowerTriangle(Matrix<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = 0; column < ColumnCount; column++)
                {
                    result.At(row, column, row > column ? At(row, column) : Zero);
                }
            }
        }

        public virtual Matrix<T> StrictlyUpperTriangle()
        {
            var result = Builder.SameAs(this);
            for (var row = 0; row < RowCount; row++)
            {
                for (var column = row + 1; column < ColumnCount; column++)
                {
                    result.At(row, column, At(row, column));
                }
            }
            return result;
        }
        public virtual void StrictlyUpperTriangle(Matrix<T> result)
        {
            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != RowCount || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            for (var row = 0; row < RowCount; row++)
            {
                for (var column = 0; column < ColumnCount; column++)
                {
                    result.At(row, column, row < column ? At(row, column) : Zero);
                }
            }
        }



        #endregion


        #region Finding

        public Tuple<int, int, T> Find(Func<T, bool> predicate, Zeros zeros = Zeros.AllowSkip)
        {
            return Storage.Find(predicate, zeros);
        }

        public Tuple<int, int, T, TOther> Find2<TOther>(Func<T, TOther, bool> predicate, Matrix<TOther> other, Zeros zeros = Zeros.AllowSkip)
            where TOther : struct, IEquatable<TOther>, IFormattable
        {
            return Storage.Find2(other.Storage, predicate, zeros);
        }

        public bool Exists(Func<T, bool> predicate, Zeros zeros = Zeros.AllowSkip)
        {
            return Storage.Find(predicate, zeros) != null;
        }

        public bool Exists2<TOther>(Func<T, TOther, bool> predicate, Matrix<TOther> other, Zeros zeros = Zeros.AllowSkip)
            where TOther : struct, IEquatable<TOther>, IFormattable
        {
            return Storage.Find2(other.Storage, predicate, zeros) != null;
        }

        /// <summary>
        /// 如果所有元素都满足谓词。 <br/>
        /// Returns true if all elements satisfy a predicate.
        /// Zero elements may be skipped on sparse data structures if allowed (default).
        /// </summary>
        public bool ForAll(Func<T, bool> predicate, Zeros zeros = Zeros.AllowSkip)
        {
            return Storage.Find(x => !predicate(x), zeros) == null;
        }

        /// <summary>
        /// Returns true if all element pairs of two matrices of the same size satisfy a predicate.
        /// Zero elements may be skipped on sparse data structures if allowed (default).
        /// </summary>
        public bool ForAll2<TOther>(Func<T, TOther, bool> predicate, Matrix<TOther> other, Zeros zeros = Zeros.AllowSkip)
            where TOther : struct, IEquatable<TOther>, IFormattable
        {
            return Storage.Find2(other.Storage, (x, y) => !predicate(x, y), zeros) == null;
        }

        #endregion


        #region COMBINATORS - Maping

        public void Map(Func<T, T> fun, Matrix<T> result, Zeros zeros = Zeros.AllowSkip)
        {
            if (ReferenceEquals(this, result))
            {
                Storage.MapInplace(fun, zeros);
            }
            else
            {
                Storage.MapTo(result.Storage, fun, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear);
            }
        }


        /// <summary>
        /// Applies a function to each value of this matrix and replaces the value with its result.
        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
        /// on the actual data storage implementation (relevant mostly for sparse matrices).
        /// </summary>
        public void MapInplace(Func<T, T> f, Zeros zeros = Zeros.AllowSkip)
        {
            Storage.MapInplace(f, zeros);
        }

        /// <summary>
        /// Applies a function to each value of this matrix and replaces the value with its result.
        /// The row and column indices of each value (zero-based) are passed as first arguments to the function.
        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
        /// on the actual data storage implementation (relevant mostly for sparse matrices).
        /// </summary>
        public void MapIndexedInplace(Func<int, int, T, T> f, Zeros zeros = Zeros.AllowSkip)
        {
            Storage.MapIndexedInplace(f, zeros);
        }

        /// <summary>
        /// Applies a function to each value of this matrix and replaces the value in the result matrix.
        /// The index of each value (zero-based) is passed as first argument to the function.
        /// If forceMapZero is not set to true, zero values may or may not be skipped depending
        /// on the actual data storage implementation (relevant mostly for sparse matrices).
        /// </summary>
        public void MapIndexed(Func<int, int, T, T> f, Matrix<T> result, Zeros zeros = Zeros.AllowSkip)
        {
            if (ReferenceEquals(this, result))
            {
                Storage.MapIndexedInplace(f, zeros);
            }
            else
            {
                Storage.MapIndexedTo(result.Storage, f, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear);
            }
        }

        public Matrix<TU> Map<TU>(Func<T, TU> f, Zeros zeros = Zeros.AllowSkip)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            var result = Matrix<TU>.Builder.SameAs(this, RowCount, ColumnCount, fullyMutable: zeros == Zeros.Include);
            Storage.MapToUnchecked(result.Storage, f, zeros, ExistingData.AssumeZeros);
            return result;
        }

        public Matrix<TU> MapIndexed<TU>(Func<int, int, T, TU> f, Zeros zeros = Zeros.AllowSkip)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            var result = Matrix<TU>.Builder.SameAs(this, RowCount, ColumnCount, fullyMutable: zeros == Zeros.Include);
            Storage.MapIndexedToUnchecked(result.Storage, f, zeros, ExistingData.AssumeZeros);
            return result;
        }

        public void MapConvert<TU>(Func<T, TU> f, Matrix<TU> result, Zeros zeros = Zeros.AllowSkip)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            Storage.MapTo(result.Storage, f, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear);
        }

        public void MapIndexedConvert<TU>(Func<int, int, T, TU> f, Matrix<TU> result, Zeros zeros = Zeros.AllowSkip)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            Storage.MapIndexedTo(result.Storage, f, zeros, zeros == Zeros.Include ? ExistingData.AssumeZeros : ExistingData.Clear);
        }



        public Matrix<T> Map2(Func<T, T, T> f, Matrix<T> other, Zeros zeros = Zeros.AllowSkip)
        {
            var result = Builder.SameAs(this);
            Storage.Map2To(result.Storage, other.Storage, f, zeros, ExistingData.AssumeZeros);
            return result;
        }

        public void Map2(Func<T, T, T> fun, Matrix<T> other, Matrix<T> result, Zeros zeros = Zeros.AllowSkip)
        {
            Storage.Map2To(result.Storage, other.Storage, fun, zeros, ExistingData.Clear);
        }





        #endregion

        #region COMBINATORS - Folding

        /// <summary>
        /// For each row, applies a function f to each element of the row, threading an accumulator argument through the computation.
        /// Returns an array with the resulting accumulator states for each row.
        /// </summary>
        public TU[] FoldByRow<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip)
        {
            var result = new TU[RowCount];
            if (!EqualityComparer<TU>.Default.Equals(state, default(TU)))
            {
                CommonParallel.For(0, result.Length, 4096, (a, b) =>
                {
                    for (int i = a; i < b; i++)
                    {
                        result[i] = state;
                    }
                });
            }
            Storage.FoldByRowUnchecked(result, f, (x, c) => x, result, zeros);
            return result;
        }

        public TU[] FoldByColumn<TU>(Func<TU, T, TU> f, TU state, Zeros zeros = Zeros.AllowSkip)
        {
            var result = new TU[ColumnCount];
            if (!EqualityComparer<TU>.Default.Equals(state, default(TU)))
            {
                CommonParallel.For(0, result.Length, 4096, (a, b) =>
                {
                    for (int i = a; i < b; i++)
                    {
                        result[i] = state;
                    }
                });
            }
            Storage.FoldByColumnUnchecked(result, f, (x, c) => x, result, zeros);
            return result;
        }

        public Vector<TU> FoldRows<TU>(Func<Vector<TU>, Vector<T>, Vector<TU>> f, Vector<TU> state)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            foreach (var vector in EnumerateRows())
            {
                state = f(state, vector);
            }
            return state;
        }

        public Vector<TU> FoldColumns<TU>(Func<Vector<TU>, Vector<T>, Vector<TU>> f, Vector<TU> state)
            where TU : struct, IEquatable<TU>, IFormattable
        {
            foreach (var vector in EnumerateColumns())
            {
                state = f(state, vector);
            }
            return state;
        }

        /// <summary>
        /// 通过在两个行向量之间应用一个函数来减少所有行向量，直到只剩下一个向量。<br/>
        /// Reduces all row vectors by applying a function between two of them, until only a single vector is left.
        /// </summary>
        public Vector<T> ReduceRows(Func<Vector<T>, Vector<T>, Vector<T>> f)
        {
            return EnumerateRows().Aggregate(f);
        }
        public Vector<T> ReduceColumns(Func<Vector<T>, Vector<T>, Vector<T>> f)
        {
            return EnumerateColumns().Aggregate(f);
        }

        /// <summary>
        /// Applies a function to update the status with each value pair of two matrices and returns the resulting status.
        /// </summary>
        public TState Fold2<TOther, TState>(Func<TState, T, TOther, TState> f, TState state, Matrix<TOther> other, Zeros zeros = Zeros.AllowSkip)
            where TOther : struct, IEquatable<TOther>, IFormattable
        {
            return Storage.Fold2(other.Storage, f, state, zeros);
        }

        #endregion


        #region Modifying

        /// <summary>
        /// Concatenates this matrix with the given matrix.
        /// </summary>
        /// <param name="right"></param>
        /// <returns></returns>
        public Matrix<T> Append(Matrix<T> right)
        {
            if (right == null)
            {
                throw new ArgumentNullException(nameof(right));
            }

            if (right.RowCount != RowCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameRowDimension);
            }

            var result = Builder.SameAs(this, right, RowCount, ColumnCount + right.ColumnCount, fullyMutable: true);
            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            right.Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, right.RowCount, 0, ColumnCount, right.ColumnCount, ExistingData.AssumeZeros);
            return result;
        }
        public void Append(Matrix<T> right, Matrix<T> result)
        {
            if (right == null)
            {
                throw new ArgumentNullException(nameof(right));
            }

            if (right.RowCount != RowCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameRowDimension);
            }

            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.ColumnCount != (ColumnCount + right.ColumnCount) || result.RowCount != RowCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameColumnDimension);
            }

            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.Clear);
            right.Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, right.RowCount, 0, ColumnCount, right.ColumnCount, ExistingData.Clear);
        }

        /// <summary>
        /// Stacks this matrix on top of the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="lower">The matrix to stack this matrix upon.</param>
        /// <returns>The combined matrix.</returns>
        public Matrix<T> Stack(Matrix<T> lower)
        {
            if (lower == null)
            {
                throw new ArgumentNullException(nameof(lower));
            }

            if (lower.ColumnCount != ColumnCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameColumnDimension, nameof(lower));
            }

            var result = Builder.SameAs(this, lower, RowCount + lower.RowCount, ColumnCount, fullyMutable: true);
            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, 0, lower.ColumnCount, ExistingData.AssumeZeros);

            return result;
        }
        public void Stack(Matrix<T> lower, Matrix<T> result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException(nameof(lower));
            }

            if (lower.ColumnCount != ColumnCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameColumnDimension, nameof(lower));
            }

            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != (RowCount + lower.RowCount) || result.ColumnCount != ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.Clear);
            lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, 0, lower.ColumnCount, ExistingData.Clear);
        }

        /// <summary>
        /// Diagonally stacks his matrix on top of the given matrix. The new matrix is a M-by-N matrix,
        /// where M = this.Rows + lower.Rows and N = this.Columns + lower.Columns.
        /// The values of off the off diagonal matrices/blocks are set to zero.
        /// </summary>
        /// <param name="lower">The lower, right matrix.</param>
        /// <returns></returns>
        public Matrix<T> DiagonalStack(Matrix<T> lower)
        {
            if (lower == null)
            {
                throw new ArgumentNullException(nameof(lower));
            }

            var result = Builder.SameAs(this, lower, RowCount + lower.RowCount, ColumnCount + lower.ColumnCount, RowCount != ColumnCount);
            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount, ExistingData.AssumeZeros);
            return result;
        }
        public void DiagonalStack(Matrix<T> lower, Matrix<T> result)
        {
            if (lower == null)
            {
                throw new ArgumentNullException(nameof(lower));
            }

            if (result == null)
            {
                throw new ArgumentNullException(nameof(result));
            }

            if (result.RowCount != RowCount + lower.RowCount || result.ColumnCount != ColumnCount + lower.ColumnCount)
            {
                throw DimensionsDontMatch<ArgumentException>(this, result, "result");
            }

            Storage.CopySubMatrixToUnchecked(result.Storage, 0, 0, RowCount, 0, 0, ColumnCount, ExistingData.Clear);
            lower.Storage.CopySubMatrixToUnchecked(result.Storage, 0, RowCount, lower.RowCount, 0, ColumnCount, lower.ColumnCount, ExistingData.Clear);
        }


        public Matrix<T> InsertColumn(int columnIndex, Vector<T> column)
        {
            if (column == null)
            {
                throw new ArgumentNullException(nameof(column));
            }

            if ((uint)columnIndex > (uint)ColumnCount)
            {
                throw new ArgumentOutOfRangeException(nameof(columnIndex));
            }

            if (column.Count != RowCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameRowDimension, nameof(column));
            }

            var result = Builder.SameAs(this, RowCount, ColumnCount + 1, fullyMutable: true);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, columnIndex, ExistingData.AssumeZeros);
            result.SetColumn(columnIndex, column);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, columnIndex, columnIndex + 1, ColumnCount - columnIndex, ExistingData.AssumeZeros);
            return result;
        }
        public Matrix<T> RemoveColumn(int columnIndex)
        {
            if ((uint)columnIndex >= (uint)ColumnCount)
            {
                throw new ArgumentOutOfRangeException(nameof(columnIndex));
            }

            var result = Builder.SameAs(this, RowCount, ColumnCount - 1, fullyMutable: true);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, 0, 0, columnIndex, ExistingData.AssumeZeros);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, RowCount, columnIndex + 1, columnIndex, ColumnCount - columnIndex - 1, ExistingData.AssumeZeros);
            return result;
        }

        public void SetColumn(int columnIndex, Vector<T> column)
        {
            if (column == null)
            {
                throw new ArgumentNullException(nameof(column));
            }

            column.Storage.CopyToColumn(Storage, columnIndex);
        }
        public void SetColumn(int columnIndex, int rowIndex, int length, Vector<T> column)
        {
            if (column == null)
            {
                throw new ArgumentNullException(nameof(column));
            }

            column.Storage.CopyToSubColumn(Storage, columnIndex, 0, rowIndex, length);
        }
        public void SetColumn(int columnIndex, T[] column)
        {
            if (column == null)
            {
                throw new ArgumentNullException(nameof(column));
            }

            new VectorStorageDense<T>(column.Length, column).CopyToColumn(Storage, columnIndex);
        }

        public Matrix<T> InsertRow(int rowIndex, Vector<T> row)
        {
            if (row == null)
            {
                throw new ArgumentNullException(nameof(row));
            }

            if ((uint)rowIndex > (uint)RowCount)
            {
                throw new ArgumentOutOfRangeException(nameof(rowIndex));
            }

            if (row.Count != ColumnCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentMatrixSameRowDimension, nameof(row));
            }

            var result = Builder.SameAs(this, RowCount + 1, ColumnCount, fullyMutable: true);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, rowIndex, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            result.SetRow(rowIndex, row);
            Storage.CopySubMatrixTo(result.Storage, rowIndex, rowIndex + 1, RowCount - rowIndex, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            return result;
        }
        public Matrix<T> RemoveRow(int rowIndex)
        {
            if ((uint)rowIndex >= (uint)RowCount)
            {
                throw new ArgumentOutOfRangeException(nameof(rowIndex));
            }

            var result = Builder.SameAs(this, RowCount - 1, ColumnCount, fullyMutable: true);
            Storage.CopySubMatrixTo(result.Storage, 0, 0, rowIndex, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            Storage.CopySubMatrixTo(result.Storage, rowIndex + 1, rowIndex, RowCount - rowIndex - 1, 0, 0, ColumnCount, ExistingData.AssumeZeros);
            return result;
        }

        public void SetRow(int rowIndex, Vector<T> row)
        {
            if (row == null)
            {
                throw new ArgumentNullException(nameof(row));
            }

            row.Storage.CopyToRow(Storage, rowIndex);
        }
        public void SetRow(int rowIndex, int columnIndex, int length, Vector<T> row)
        {
            if (row == null)
            {
                throw new ArgumentNullException(nameof(row));
            }

            row.Storage.CopyToSubRow(Storage, rowIndex, 0, columnIndex, length);
        }
        public void SetRow(int rowIndex, T[] row)
        {
            if (row == null)
            {
                throw new ArgumentNullException(nameof(row));
            }

            new VectorStorageDense<T>(row.Length, row).CopyToRow(Storage, rowIndex);
        }

        public virtual void SetDiagonal(Vector<T> source)
        {
            if (source == null)
            {
                throw new ArgumentNullException(nameof(source));
            }

            var min = System.Math.Min(RowCount, ColumnCount);

            if (source.Count != min)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentVectorsSameLength, nameof(source));
            }

            for (var i = 0; i < min; i++)
            {
                At(i, i, source.At(i));
            }
        }
        public virtual void SetDiagonal(T[] source)
        {
            if (source == null)
            {
                throw new ArgumentNullException(nameof(source));
            }

            var min = System.Math.Min(RowCount, ColumnCount);

            if (source.Length != min)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentArraysSameLength, nameof(source));
            }

            for (var i = 0; i < min; i++)
            {
                At(i, i, source[i]);
            }
        }


        public void SetSubMatrix(int rowIndex, int columnIndex, Matrix<T> subMatrix)
        {
            subMatrix.Storage.CopySubMatrixTo(Storage, 0, rowIndex, subMatrix.RowCount, 0, columnIndex, subMatrix.ColumnCount);
        }
        public void SetSubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount, Matrix<T> subMatrix)
        {
            subMatrix.Storage.CopySubMatrixTo(Storage, 0, rowIndex, rowCount, 0, columnIndex, columnCount);
        }
        public void SetSubMatrix(int rowIndex, int sorceRowIndex, int rowCount, int columnIndex, int sourceColumnIndex, int columnCount, Matrix<T> subMatrix)
        {
            subMatrix.Storage.CopySubMatrixTo(Storage, sorceRowIndex, rowIndex, rowCount, sourceColumnIndex, columnIndex, columnCount);
        }

        #endregion


        #region Transpose

        public Matrix<T> Transpose()
        {
            var result = Builder.SameAs(this, ColumnCount, RowCount);
            Storage.TransposeToUnchecked(result.Storage, ExistingData.AssumeZeros);
            return result;
        }
        public void Transpose(Matrix<T> result)
        {
            Storage.TransposeTo(result.Storage, ExistingData.Clear);
        }

        // 共轭转置
        public virtual Matrix<T> ConjugateTranspose() => throw new Exception();
        public virtual void ConjugateTranspose(Matrix<T> result) => throw new Exception();

        /// <summary>
        /// 根据排列来计算矩阵的行 <br/>
        /// Permute(序列改变) the rows of a matrix according to a permutation.
        /// </summary>
        /// <param name="p">The row permutation to apply to this matrix.</param>
        public virtual void PermuteRows(Permutation p)
        {
            if (p.Dimension != RowCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentArraysSameLength, nameof(p));
            }

            // Get a sequence of inversions from the permutation.
            var inv = p.ToInversions();

            for (var i = 0; i < inv.Length; i++)
            {
                if (inv[i] != i)
                {
                    var q = inv[i];
                    for (var j = 0; j < ColumnCount; j++)
                    {
                        var temp = At(q, j);
                        At(q, j, At(i, j));
                        At(i, j, temp);
                    }
                }
            }
        }
        public virtual void PermuteColumns(Permutation p)
        {
            if (p.Dimension != ColumnCount)
            {
                throw new ArgumentException(Resources.Numerics.ArgumentArraysSameLength, nameof(p));
            }

            // Get a sequence of inversions from the permutation.
            var inv = p.ToInversions();

            for (var i = 0; i < inv.Length; i++)
            {
                if (inv[i] != i)
                {
                    var q = inv[i];
                    for (var j = 0; j < RowCount; j++)
                    {
                        var temp = At(j, q);
                        At(j, q, At(j, i));
                        At(j, i, temp);
                    }
                }
            }
        }


        #endregion


        //}}@@@
    }


}

